clear; clc; close all;

%% 固定参数（除 b0, b1 外）
a0 = -2;
a1 = 0.83;
a2 = 0.14;
b2 = 1;
d1 = 0;
d2 = 0;

% 扫描参数 b0 和 b1 的取值范围
b0_values = linspace(-5, 5, 101);
b1_values = linspace(-1, 1, 101); 
nb0 = length(b0_values);
nb1 = length(b1_values);

% 迭代设置
N_total = 2000;    % 总迭代步数
N_trans = 1000;     % 舍弃瞬态步数
N_show  = N_total - N_trans;  % 有效记录步数

% 初始条件
x0 = 0.1;
y0 = 0.1;
z0 = 0.1;

%% 预分配存储 C0 复杂度矩阵，尺寸为 (nb1 x nb0)
% 约定外层循环为 b0, 内层为 b1
C0 = zeros(nb0, nb1);

%% 双重循环，遍历每个 (b0, b1) 组合
for i = 1:nb1
    for j = 1:nb0
        b1 = b1_values(i);
        b0 = b0_values(j);
        
        % 初始化状态，采用初始条件
        x = x0; y = y0; z = z0;
        
        % 舍弃瞬态迭代
        for n = 1:N_trans
            [x, y, z] = vl_dmsin(a0, a1, a2, b0, b1, b2, d1, d2, x, y, z);
        end
        
        % 记录后续 N_show 步的 x 值到 data 数组（用 x 作为复杂度计算依据）
        data = zeros(1, N_show);
        for n = 1:N_show
            [x, y, z] = vl_dmsin(a0, a1, a2, b0, b1, b2, d1, d2, x, y, z);
            data(n) = x;
        end
        
        C0(i,j) = C0_Complexity(data, 1);
    end
    fprintf('%.1f%%\n',100 * i / nb1);
end

%% 生成绘图所需的网格（注意行列对应关系，外层 b1, 内层 b0）
[B0, B1] = meshgrid(b0_values, b1_values);

%% 绘图
figure;
surf(B1, B0, C0);
shading interp;
xlabel('b_1','FontSize',12);
ylabel('b_0','FontSize',12);
zlabel('C_0 复杂度','FontSize',12);
title('Fig.19(b)： C_0 复杂度, b_0, b_1','FontSize',14);
colorbar;


function c0_norm = C0_Complexity(data, energyThreshold)
% C0_Complexity_v2 计算归一化的 C0 复杂度（基于 FFT 能量分解方法）
%
%   c0_norm = C0_Complexity_v2(data, energyThreshold)
%
% 输入参数:
%   data           - 时间序列（向量），例如来自离散映射的一段数据
%   energyThreshold- 能量累积阈值（介于 0 到 1 之间），用于选择主要频率成分；
%                    默认取 0.5，即保留能量累积达到 50% 的频率成分
%
% 输出参数:
%   c0_norm        - 归一化的 C0 复杂度，数值范围 [0, 1]
%
% 算法说明：
%   1. 对数据去均值后做 FFT，获得能量谱；
%   2. 按能量从大到小对频率成分进行排序，找出最少的频率数 k，使得它们的能量累计
%      占总能量的至少 energyThreshold；
%   3. 将剩余的频率成分置零，再做逆 FFT 得到“规则信号”；
%   4. 定义残差信号为原始信号与规则信号的差，其能量占总能量的比例作为 C0，
%      如果信号越难以用少数频率描述，则残差能量越大，C0 越接近 1；反之越接近 0。
%
% 注意：本函数中，归一化时直接用残差能量除以原信号能量，确保输出值在 [0,1] 范围内。

    if nargin < 2
        energyThreshold = 0.5; % 默认阈值 0.5，可以根据需要修改
    end

    data = data(:)';  % 转为行向量
    N = length(data);
    
    % 去中心化数据
    data_centered = data - mean(data);
    
    % FFT 计算
    Y = fft(data_centered);
    P = abs(Y).^2;  % 能量谱
    totalEnergy = sum(P);
    
    % 对能量谱排序，获取排序索引
    [sortedEnergy, sortIdx] = sort(P, 'descend');
    cumEnergy = cumsum(sortedEnergy);
    
    % 找出最小 k，使得累计能量达到 energyThreshold * totalEnergy
    k = find(cumEnergy >= energyThreshold * totalEnergy, 1, 'first');
    
    % 构造频谱掩模：保留排序前 k 个频率成分，其余置零
    mask = zeros(size(Y));
    mask(sortIdx(1:k)) = 1;
    
    % 得到主要频率信号
    Y_main = Y .* mask;
    ruleSignal = real(ifft(Y_main)) + mean(data);  % 加回均值
    
    % 计算残差能量
    residual = data - ruleSignal;
    resEnergy = sum(residual.^2);
    
    % 定义 C0 为残差能量占总能量（也可根据需要使用 data_centered 的总能量）
    c0 = resEnergy / sum(data.^2);
    
    % 为确保数值在 [0,1] 内，可适当裁剪
    c0_norm = min(max(c0, 0), 1);
end

